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ABSTRACT 
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The  numerical  construction  of  a  smooth  surface  with  prescribed  weighted 
integrals  over  a  domain  of  interest,  is  investigated.  This  construction  is 
mostly  relevant  to  the  estimation  of  a  smooth  density  function  over  qeoqraphical 
regions,  from  data  aggregated  over  several  subregions.  By  analoqy  to  the  defini¬ 
tion  of  the  univariate  histospline  the  smooth  surface  is  defined  as  the  solution 
to  a  certain  constrained  minimization  problem.  The  application  of  finite  element 
methods  to  the  numerical  solution  of  this  minimization  problem  is  studied.  It 
is  shown  that  any  finite  element  procedure,  convergent  for  a  related  boundary 


value  problem  can  be  used  to  construct  a  sequence  of  finite  element  approximations 
converging  to  the  smooth  surface  which  solves  the  constrained  minimization  problem. 

For  the  case  of  smoothness  requirement  of  lowest  order,  a  specific  finite 
element  method  is  considered,  and  its  converqence  as  the  mesh  size  decreases  is 
demonstrated  numerically  for  a  particular  example  of  ^volume  matchinqf'. 
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SIGNIFICANCE  AND  EXPLANATION 


We  consider  a  numerical  method  for  the  construction  of  a  smooth  donsily 
function  from  available  data  in  aggregated  form.  For  example,  suppose  th.it  the 
population  census  is  given  by  bureaucratic  region  (say,  states')  and  it  is  desir.  ,1 
to  construct  numerically  a  smooth  function  f(x,y)  intended  to  be  an  estimate  o! 
the  population  density  at  location  (x,vj.  In  order  to  select  from  the  infinit.lv 
many  ways  in  which  this  could  be  done  a  particular  one,  we  require  that  f  bo  the 
"smoothest"  function  matching  the  prescribed  aggregated  data.  Our  measure  of 
roughness  to  be  minimized  by  f  is  the  integral  over  the  region  of  interest  of  a 
quadratic  form  in  all  the  derivatives  of  f  of  a  certain  order.  This  order  is  a 
free  parameter  which  can  be  chosen  according  to  the  required  degree  of  smoothness. 
By  using  finite  element  techniques  we  reduce  the  computation  of  the  minimal  f  to 
that  of  solving  a  finite  dimensional  constrained  minimization  problem  of  a  particu¬ 
lar  structure.  We  show  that  any  finite  element  scheme  which  produces  good  approxi¬ 
mations  to  the  solution  of  a  related  elliptic  boundary  value  problem  can  be  used 
in  order  to  produce  good  approximations  to  the  required  smooth  surface.  This 
method  is  discussed  in  detail  for  the  particular  case  of  "volume  matching",  under 
the  requirement  of  minimal  integral  of  the  sum  of  squares  of  the  first  partial 
derivatives . 
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NUMERICAL  CONSTRUCTION  OF  SMOOTH  SURFACES  FROM  AGGREGATED  DATA 


*  ** 

Nira  Dyn  and  Warren  Ferguson 


1.  Introduction 

This  work  is  motivated  by  the  need  for  a  numerical  procedure  for  the  construction  of  a 
smooth  surface,  describing  a  certain  geographically  varying  quantity  over  a  finite  geographi¬ 
cal  region,  given  the  integrals  of  the  quantity  over  several  disjoint  sub-regions.  One  of 
these  problems  is  that  of  estimating  the  density  of  a  population  over  an  area  as  a  smooth 
function  of  the  geographical  coordinates,  given  the  population  census  according  to  a  certain 
bureaucratic  subdivision  of  the  area  [11]  ,  [S] .  In  this  context  the  additional  constraint  of 
positivity  of  the  surface  is  in  place. 

A  method  for  estimating  a  multivariate  smooth  function  from  aggregated  data  is  presented 
and  analyzed  in  [5] .  The  function  is  chosen  by  minimizing  a  region-dependent  roughness 
criterion  subject  to  the  given  aggregated  data.  For  the  bivariate  case  and  for  homogeneous 
roughness  criteria  the  estimating  surface  is  taken  as  the  solution  of  the  following  minimiza¬ 
tion  problem: 

Problem  I:  Find  u* 

m  f  ,m  \  2 

(1.1)  minimizing  J  (u)  =  f  V  (m)  ■ — -  U  .  t  dxdy 

m  J  t  i  '  i  m-i I  ’ 

si  1=0  \dx  3y  / 

among  all  function  satisfying 

(1.2)  L.u  :  /  u<£ .  =  s.  ,  i  =  l,...,N 

i5 

2  2 

Mere  is  a  bounded  domain  in  R  ,  $  ,  L  (>:)  i=l,...,N  and 
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,k 

Hm(i7)  =  { u  |  - r-  -  £  L  (fi)  ,  0<iik,  0  <  k  <  m} 

ax^y14'1 


Since  J  (u)  vanishes  on  O  -  the  (mt*>  dimensional  space  of  all  polynomials  of  total 
m  ~m  2 

degree  <  m,  this  method  of  estimating  a  surface  reproduces  any  surface  which  is  a  poly¬ 
nomial  in  Q^,  whenever  there  is  only  one  polynomial  in  0^  satisfying  the  constraints  (1.2). 
Thus  the  degree  of  smoothness  in  m,  which  is  a  free  parameter,  can  be  chosen  according  to  the 
required  smoothness  properties  of  the  surface,  but  with  the  obvious  limitation  that  0  does 


not  contain  c.  nontrivial  polynomial  satisfying  (1.2)  with  s  =‘*-=s  =  0.  In  particular  this 

1  N 

implies  that  l"1*1)  IN.  This  approach  is  similar  to  that  of  [4),  [9],  where  interpolating 
surfaces  are  constructed  by  minimizing  region  independent  roughness  criteria  of  the  form 


m  /  ,m  Y 
l  <”>  I 

i=0  1  \3x18ym  7 


The  solution  to  Problem  I  is  characterized  in  [5] ,  and  is  shown  to  be  related  to  a  certain 
elliptic  boundary  value  problem.  This  solution  can  be  regarded  as  a  generalization  of  the 
concept  of  univariate  histosplines  [21 .  The  "volume  matching"  problem  in  a  tensor-product 
situation  is  studied  in  [10J ,  where  the  solution  is  shown  to  be  a  tensor-product  of  uni¬ 
variate  histosplines  and  where  a  computational  algorithm  is  presented. 

Tn  the  present  work  we  investigate  the  applicability  of  finite  element  methods  developed 
for  the  solution  of  elliptic  boundary  value  problems,  to  the  construction  of  approximations  to 
tne  solution  of  Problem  I.  We  discretize  Problem  I  by  minimizing  (1.1)  subject  to  (1.2)  among 
all  functions  in  a  finite  dimensional  subspace  of  II  (  j  spanned  by  "finite  elements".  It 
is  shown  that  any  finite  element  scheme,  convergent  for  a  related  elliptic  boundary  value 
problem,  can  be  used  to  construct  a  sequence  of  finite  clement  approximations  converging  to 
the  solution  of  Problem  I. 

In  case  the  surface  is  required  to  be  nonnegative  we  are  led  to: 

Problem  II:  find  u  ■  nm(  )  minimizing  (1.1)  among  al)  nonnegative  functions  satisf/mu 

(1.2)  . 

for  the  case  m  =  1  the  discretized  version  of  Problem  II  becomes  a  quadratic  pro¬ 
gramming  problem  when  the  discretization  is  made  by  nonnegative  finite  elements,  and  a 
characterization  of  the  discretization  solution  is  given. 


-2- 


The  numerical  computation  of  the  solution  to  the  discretized  Problem  I  involves  a  solution 
of  a  large  scale  linear  system  of  a  special  structure,  with  a  main  part  resemblinq  the  linear 
system  characteristic  to  the  finite  element  solution  of  the  related  boundary  value  problem. 
Iterative  schemes  for  the  computation  of  the  solution  of  such  large  systems  are  analyzed  m 
[6]  . 

/  2  2  \  2 

Problem  II  is  formulated  for  m  =  1  and  for  J( u)  -  f  I-1— y  ♦  — rj  ^xc*y  in  the  context 

V*  3y  / 

of  "volume  matching"  in  (11)  where  an  iterative  procedure  is  presented  for  the  numerical 
computation  of  an  approximate  solution.  The  procedure  is  tested  on  several  examples,  and  the 
iterations  converge.  Yet  no  proof  of  convergence  is  given.  This  procedure  with  the  additional 
condition  u  -  0  on  the  boundary  T  of  and  without  the  steps  that  impose  nonneaat i vi t y , 

is  in  fact  one  of  the  converging  iterative  schemes  of  161  (see  also  Section  5)  .  The  inj  1  orien¬ 
tation  of  the  nonnegativity  in  Ill)  seems  to  be  incorrect,  in  view  of  the  character izat ion  of 
the  solution  to  the  discretized  version  of  Problem  II  (Section  \) ,  and  not  alwavs  possible  as 
remarked  in  (11). 

Iterative  schemes  for  the  computation  of  solutions  to  the  discretized  version  of  Problem 
II  for  m  =  1,  and  the  convergence  of  these  solutions  to  the  solution  of  the  cont  muon*, 
problem,  are  yet  under  investigation. 

In  Sect  ion  the  results  in  |'d  concerning  the  *-.olut  ion  of  Ptoblem  1  ait*  leviewed.  In 

Section  *  we  discretize  Problem  !  and  charactei  i/.e  the  solution  of  the  di-vietiz.d  ptoblem 

for  m  l»  and  for  m  ~  1  with  the  nonneg.it  i  v  1 1  v  constraints.  The  convergence  of  the  finite 
element  solutions  l  s  dealt  with  jn  Section  -1,  wl.i  h  \  com)  ut  at  iona  1  method  im  t  ho  numeti  al 
solution  of  t  Ik  dl  set  rt  i  zed  Jiobliin  n  1,  uhd  it  mu  (orient  at  i « *i.  on  a  .-ompitet  t  u  ts« 

"volume  matching"  ptoblem,  ace  d «  ic  •  n.  Stvt  t«>n  •>.  The  convet  uetice  of  this  s.-hen*  a  t  s> 

mesh  •.  i  >!•>  i  t  us.  i  •.  d«  *rt*  -n  t  »  a  t  *  d  nun*  i  i  ca  1  1  •  !  <  o  a  i  .n  t  i.  ul.ti  *  xapi}  1  e  ot  "  Vo  1  urn*  mat  .  sine". 


Characterization  of  the  solution  of  Problem  I. 


The  analysis  of  Problem  I  ir.  fr,l  is  made  utvU  r  a  i  .»*  \  .  ;  ;r-;  * 

ir\+  1 

of  (l.i)  with  the  subspace  0  of  dimension  M  -  i  )  ,  o >:.*:*-♦  i: 

*-m 

m  x , y  of  total  decree  <  m. 

Assumption  M. 

The  only  polynomial  o  .  o  which  satisfies  f  a  $ .  -  ,  i  i,...,:, 

m  *  i 

equivalently 


,  r  .  m  #  n 

rank'  !  q  t  .  ,  .  ,  *  M,  wit*. 

J  )  i*l,3»l 


The  characterization  of  the  solution  or  I  i  oi  l«  n  !  ; 
form  associated  with  the  functional  J  of  <  1  . 1 '  : 


'  *•*  m  •  .1  •  v 

A  { u / V' )  =  t  '  { i  -  -  -  -:x.: 

n  '  •  i  l  r.-i  :  r- 1 

-x  *•  •>: 


It  is  shown  ir.  I )  that  the  solution  o:  lie:  i. 


sat  i  M  vinn  : 


A  ( u * » v )  -  -  i  *  .  :  )  v  :  oi  .i  1 1 

m  *  i*l  1  * 


L  u*  :  u*: 

i  •  i 


constants.  Since  A  tv,q)  -  !oi  i'.  . 


.  i  .  .  :.ui  a  1 1  i  i  c.»?  1 1  *r.  * 


i  «  ’*  l  •.  t  :.t  j  : .  l  :  \  *  >  l  •  n>  *m  i .» 1 


12.7) 


i«l , . . .  ,M 


L.q*  5  fq*4>.  =  s.  , 

1  \  1  1 

and  where  ,  1  i  <  N-M,  is  the  unique  solution  in  Hm(.q  of  the  following  boundary  vui 
problem,  formulated  variationally  as: 

(2.8)  A  (5.,v)  =  f  $.v  for  all  v  t  Hm ( .1) 

m  r  i 

(2.9)  |  r  i<t>  =  o,  j=l,...,M  . 


In  (2.8) 

(2.10) 


M 

d .  =  6 .  -  \  \  ,  g . 

V1  i+M  11  1 

■,  =  1  J 


where  are  constants  determined  by  the  condition: 

(2.11)  f  qg .  =  0  for  all  q  t  Q 

SI  1 

In  case  11  is  a  smooth  domain  the  boundary  value  problem  (2.8),  (2.9)  can  he  reform:.:  In 
as  (51  : 

( 32  32\  - 

!  — r  +  — -  'f, .  =  <> .  in 

\,2  2  .  1  l 

\3x  3y  / 

(2.12)  S^£^=0onr  m  <  j  2m- ) 

Mi  =  I  4^,  =  3  =  1 . M  • 

In  (2.12)  T  is  the  boundary  of  i',  and  6  .  are  differential  o:orat>u:  '!  : 

m  2m- 1  1 

m,...»2m-l,  such  that  the  generalized  Green  Formula  holds: 


m  r 

.i-  * 

*  (-1) 

J 

V  — 

+  — 

i.' 

X 

m- 1 

+  y 

/ 

A 

u  — 

j’O 

’r 

2m- 1  - 

'3  in 

with  — r  the  jth  normal  derivative  of  v. 

an3 


3 .  Discretization  of  Problem  I  and  Problem  II. 

Let  be  the  span  of  n  linearly  independent  functions  v  , satisfying  the 

following  assumption: 

Assumption  3.1. 

V  contains  the  space  0  ,  and 


the  linear  functionals  L,,...,L„  in  (1.2)  are  linearly  independent  over  V  . 

IN  n 

The  first  requirement  in  Assumption  3.1  can  be  met  with  v^,...,v  piecewise  polynomials 

in  x,y  of  total  degree  <  m,  with  local  supports,  such  that  \r  t  cm  *(2),  i=l,...,n.  The 

second  requirement  can  be  guaranteed  by  taking  n  large  enough  and  the  supports  of  v^,...,v 

small  enough,  with  their  union  containing  P. .  In  the  "volume  matching"  problem,  where 

(3.1)  L^u  =  J  u,  i=  1 , . . . , N 


for  s) . il  a  partition  of  .0,  the  second  part  of  Assumption  3.1  is  satisfied,  if  to  each 

1  N 

there  corresponds  a  nonnegative  function  in  whose  support  is  contained  in  .  j,.  Both 

requirements  can  be  met  if  v^,...,v  (n'N)  are  tensor  products  of  mth  degree  univariate 

B-splines  [10],  although  such  a  choice  is  not  efficient  computationally  for  n  -  1  . 

The  discretized  version  of  Problem  I  with  respect  to  V  is  13): 


Problem  I:  Find  v  i  V 
- D  n 


minimizing  J  (v) 


among  all  functions  satisfying  (1.2). 

Theorem  3.1:  There  exists  a  uniuue  v  ■  V  which  solves  Problem  I  .  This  solution  is 
— — — -  •  n  D 

the  uniuue  element  in  V  which,  satisfies: 


(3.2) 


A  ( v  ,  v  .  )  - 


i  —  l , 


i-1 


(3.3) 


f  v: 


l  =  1  ,  .  .  .  ,N 


whc-r*’  ^  i  jro  cnnsLints  s.it  i  nf  vi  jm 

3  1  =  1 
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(3.4) 


£  yn  /  q<f  ■  =  0  for  all  q  l  Q 


j=l  J  S3  3 


Proof:  Any  solution  of  Problem  1^,  v  =  V  ^ivi'  mus1:  satisfy  the  necessary  conciit 

i=l 

N 

(3.5)  — —  (J  (v)  +  y  \n[J  v$ .  -  s.)}  =  0  i=l,...,n 

dt> .  m  .  ,  3  -  1  1 

i  3  =  l  i2J 

since  Problem  IQ  is  an  n-dimensional  constrained  minimization  problem.  In  (3.5)  1 -j  - 

are  the  Lagrange  multipliers  corresponding  to  the  constraints  (3.3).  Dif ferentiating  (3 
we  obtain  conditions  (3.2).  In  view  of  the  first  part  of  Assumption  3.1,  (3.4)  follows 
directly  from  (3.2)  . 

In  order  to  demonstrate  the  existence  of  a  unique  element  v  in  satisfying  th 

necessary  conditions  (3.2)  and  (3.3),  it  is  sufficient  to  show  that  the  system 


A  e\,V\  /?' 

E-  oJ[yJ  [s, 


admits  a  unique  solution,  where  b  =  (b,  , . .  .  ,b  )  1  ,  y  =  ,v)  ’  ,  s  =  (s,  ,  .  .  .  ,s  ,)  ’  a 

"I  h  ~  1  N  ~  1  .\ 

the  matrices  A  and  E  are  defined  as 

(3.7)  A  =  {A  (v. ,v .) }n  ,  ,n  ,  E  =  {/  v.$.}n  ,N  . 

nxn  m  i  3  i=l, 3=1  n*N  \  13  1=1,3=! 


Thus  to  complete  the  proof  we  show  that  the  system 


A  E  \  /  x 


E*  0/\y 


n  N 

with  x  '  R  ,  y  <  R  ,  admits  only  the  trivial  solution.  Now  the  matrix  A  is  nonr.c-at  .  \ 


.'finite  since  by  (3.7)  and  (2.2)  for  any  x  c  R 


T 

x  Ax  '  A 


( ?  S  \ 

'  X , v . ,  '  XV-  '  0  Wl 

1  .'ll  1  1  ]  )!- 


m  1  .  i  1  i 

\i  =  l  ]=1 


th  eciualitv 


iff  y  x . v .  ,  0 

. 11  m 
i=l 


Rewriting  (3.8)  we  obtain 


(3.9)  Ax  +  Ey  =  0,  E'x  =  0 

n 

and  therefore  x'Ax  =  -x'Ey  =  0,  which  together  with  (3.8)  implies  that  ]_  x.v.  *  Q  .  Or. 


the  other  hand  the  equation  E'x  =0  of  (3.9),  in  view  of  (3.7),  becomes 


/  n  \ 

J  l  VjK  =  °>  i=1' 

P.  \j«l  1  7 


Since  )  x.v.  c  0  ,  (3.10)  is  consistent  with  Assumption  2.1  only  if  x  =  0.  Hence  fcv  (3.9) 

j-1  3  3  m  N  *  ' 

Ey  =  0  or  equivalently  J  v.  £  yd.  =  0,  i=l,...,n  which  is  consistent  with  Assumption  3.1 

n  1  j=i  3  3 

only  if  y  =  0. 

When  the  functions  vi'""''vn  are  local  support,  as  is  the  case  ir.  the  finite  element 

method,  the  nonnegative  definite  matrix  A  in  (3.6)  (and  in  (3.15))  is  sparse,  and  of  special 

structure.  These  properties  of  A  allow  one  to  efficiently  solve  the  system  (3.6)  using  or.e 

of  the  iterative  methods  analysed  in  (6).  (One  such  method  is  presented  in  Section  5). 

In  analogy  to  the  characterization  (2.6)— (2.11)  of  the  solution  of  Problem  I  we  have: 

n 

Theorem  3.2.  Let  b,y  be  the  solution  of  (3.6).  Then  v  =  b.v.  has  the  unique  represor.t- 
N-M  ~  ~  1=1  1  1 

ation  v  =  £  c.ry  +  q* ,  with  q*  e  Q  satisfying  /  at  .  «  s .  ,  1=1,...,.",  and  where  for 

i=l  .2 

1  <.  i  <  N-M,  t  Vn  is  determined  uniquely  by 

(3-u>  wv  =  I  Vr  j=i . n  • 


f  .  =  o,  j=i, . . . ,m  , 


with  defined  by  (2.10). 

Proof:  First  we  observe  that  (3.11)  and  (3.12)  admit  a  unique  solution.  In  fact  the  I;:.-. 

system  corresponding  to  (3.11)  is  singular,  but  by  (3.8)  and  (2.10)  the  right-hand  Sid.  ,-f 

system  is  in  the  span  of  the  columns  of  the  matrix  of  the  system  -A.  Hence  (3.11)  admit 

solution,  and  all  solutions  differ  by  a  polynomial  in  O  .  Conditions  (3.12)  determine  .= 

m 

unique  element  from  this  set  of  solutions,  in  view  of  (2.1). 
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The  functions  n, f-..*n  are  linearly  independent  in  view  of  (3.11)#  (2.10)  and  the 
1  N-M 

linear  independence  of  L.  , . * . , L  over  V  .  Therefore  the  (H-M)  *  (IJ-M)  matrix  with  entries 
*  IN  n 

r  .  N-M 

{*  :'i .  y .  '•*.  .  ,  is  non-sincrular ,  since  bv  (3.12),  (2.10)  and  (3.11) 

lvj+M  i,j*l 


ir  j+M 


i  3 


(n. 
*  3 


n . ) 

i 


and  the  bilinear  form  A  (•,*)  is  an  inner-product  in  the  subspace  (u|u  *:  H*”  (!'.), 
m 

f  u>)^  -  0,  .  Hence  aiven  s,  there  exists  a  unique 

c  *  1  such  that 

1  N-M 

/N~M  \ 


(3.13) 


j“l 


3  3, 


V  . 

1 


s i »  i  =  M+l , . . . ,N  , 


N  n 

and  the  function  v  =  V  c.n.  +  q*  =  V  b  v.  satisfies  (3.2)  ana  (3.3),  in  view  of  (3.11), 

.‘-,11  . i  1 

1=1  i=l 

(3.12)  and  (3.13).  This  completes  the  proof  of  the  theorem. 

Before  proceeding  to  the  analysis  of  the  convergence  of  the  solutions  of  a  sequence  of 
discretized  problems  to  the  solution  of  the  continuous  problem,  we  consider  the  discretized 
version  of  Problem  II.  The  discretization  is  done  by  considering  Problem  II  in  the  subspace 

.  ,  d.v.  >  0  in  .“. ,  minimizes 

r.  in  -1=1  11  — 

J^lu)  of  (1.1)  among  all  positive  functions  in  satisfying  (1.2).  In  case  the  condition 

n 

^  d.v.  >0  in  is  equivalent  to  the  condition  d  >  0,  the  discretized  version  of  Problem 

i=l  1  1 

II  becomes  a  quadratic  programming  problem  in  terms  of  A,  E,  s,  d,  namely: 

Problem  II  Find  a  <  Rn 


(3.14)  minimizing  d’Ad 

among  all  vec.  irs  satisfying  E'd  =  s,  d  >_  0 

The  solution  of  (3.14)  is  characterized  as  the  unique  solution  to  the  problem  [8]: 

Ad-EY  0 
d'(Ad-EY)  =  0 

(3.15) 

E'd  =  s 
d  >  0 
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with  >  =  (•> . >  )  1 


N  constants  determined  uniquely  by  (3.15). 


For  the  case  m»l  the  functions  in  H  (:?)  are  not  necessarily  continuous  (Nm(£)  r  C(..) 
for  m  2l  2) ,  and  the  positivity  must  be  interpreted  not  pointwise  but  in  the  following  sense: 
u  t  H^(£)  is  “positive*'  in  il,  if  it  is  the  limit  of  a  sequence  of  positive  functions  in 
C  (i'.)  in  the  Sobolev  norm 

|i  2  r  2  3u  .  3u  1  J  - 

iM  i  a  J  lu  +  {—)  +  (r-)  1  dxdy  . 

h^o)  ax  3v 


In  this  case  <-  i)  one  can  take  v^ . to  be  piecewise  linear  with  local  supports,  such 


that  v.  >  0  in  si,  and  that  for  a  reqular  mesh  of  points  (x 


:,  v  (x  )  =  '  , 

il  il 


i.j=l..--.n.  Such  a  choice  is  furnished,  for  example,  by  supports  of  the  form: 


Correspond ino  to  each  support  of  the  above  form,  centered  at  x^,  the  finite  element  vv  is 
a  linear  function  within  its  support  which  vanishes  on  the  boundaries  of  the  support  and  satis¬ 
fies  v.(x  )  *  1,  v.(x)  •  v.  For  such  a  choice  of  finite  elements 

li  l 


(3. 16) 


n 

)  d  .  v  .  ■  0  i  n 

‘•li 


d  =  (d,  , .  .  .  ,d  >  1  •  0 
1  n  -  - 


and  the  characterization  li.lfd  is  valid.  This  character izat ion  is  the  key  to  the  dovelo|-ment 
of  an  efficient  alqorithm  for  the  solution  of  (1.1-1)  for  a  larqe  and  sparse  matrix  A,  win  eh 
takes  into  account  the  special  structure  of  A.  We  intend  to  invest  mate  this  algorithm  else¬ 
where  . 
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4 .  Convergence  of  the  discretized  solutions. 

The  characterization  (2.3),  (2.4)  of  the  solution  to  Problem  I  differs  from  the  cr.aractt 
ization  of  weak  solutions  to  elliptic  boundary  value  problems  11),  (31,  by  the  fact  that  tiv 
right-hand  side  of  (2.3)  is  not  qiven,  but  instead  its  structure  up  to  !,’  constant?  is  know: 
and  there  is  additional  information  on  the  solution  in  integral  form  qiven  by  (2.4).  Never¬ 
theless  in  the  analysis  of  the  convergence  of  a  finite  element  scheme  to  the  solution  of 
Problem  I,  the  convergence  of  that  finite  element  scheme  to  the  solution  of  the  following 
boundary  value  problem  is  relevant: 

(4.1)  Am(u,v)  =  /  fv  for  all  v  e  Hm(P.),  f  c  L2(fl),  /  fq  =  0,  q  e  Q  , 

j)  n 

L^u  =  0,  i=l,  . . .  ,M 

Let  denote  a  space  of  finite  elements  such  that  the  area  of  the  support  of  each 

2  2 

element  is  bounded  from  below  and  above  by  ah  and  8h  respectively  (0  <  a  •-  ■)  .  The 

finite  elements  approximation  to  the  solution  of  (4.1)  is  u.  e  V,  satisfying  (1),  13): 

h  n 


(4.2) 


A  (u.,v)  =  j  fv,  v  e  V 
m  n  h 

L.u,  »  0  i=l , .  . .  ,M  . 
i  h 


Theorem  4.1.  Let  u,u^  be  determined  by  (4.1)  and  (4.2)  respectively,  and  let  v._ 
the  finite  elements  approximation  to  the  solution  of  Problem  I,  u*.  If 


(4.3)  ||u-u  ||  B(u)hV,  v  >  0 

where  B(u)  is  a  constant  depending  on  u,  and  |]  •  j|  is  a  norm  with  the  proyertv 


(4.4) 


-V 


v  «  Hml ') 


L  (>!) 


then  for  h  N  0  small  enough 


(4.5) 


!vh-u*| 


Gh 


with  G  =  G(n,-i'- 


1 


N  1 


.  s  ) 
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Proof:  By  the  result  cited  in  Section  2  in  (2.2)  -  (2.6) 


Therefore  V  <  (U-M)G  h  ,  and  for  h  small  t'noucih  such  that 
*  ~  1 

yields 


-1 


( 4  .  UM 


,  -1 


14.14^ 


G.  h 


with  <J,  *  2;'  L  V 


1  ’ 


To  complete  the  [  roof  of  the  theorem  observe  that  by'  <4.0  and  (4.H) 

N-M  N-M  N-M  N-M  M-M 

v  -u*  «  J  eh  -  l  c*i.  «  '  c*(n  ^  <c  -c*)C  )  +  '>  <c 

h  ,ii  “  ii  -  ill  ,ii  ii  ,iii 

i=l  i=l  i=l  i=l  1=1 


which  in  view  of  the  bounds  (4. it)  and  (4.14)  becomes: 

N-M  N-M  N-M 

.s  V  ,-u*',  •  '  i  C  *  ,  B  (:,  )  h '  *  v  U  h'B(r.  )h‘  ♦  v  U  h'  '  V.il  ’  , 

1=1  111  i*l  1  1  1=1  ‘  1  ~ 

with  <J  depending  on  the  bilinear  form  A  (•,•),  the  functions  ;  ,  and  the  data 

m  1  M 

si . V 

Thus  we  have  reduced  the  convergence  problem  related  to  the  solution  of  Problem  I  by 
finite  elements  schemes,  to  the  convenience  of  those  schemes  for  elliptic  boundary  value 
problems  of  the  form  (4.1)  -  the  Neumann  problems. 
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5.  Computation  of  a  volume  matching  surface  for  m»l. 


In  this  section  we  specialize  to  the  volume  matchinq  problem  with  m-1.  Problem  I  then 
becomes : 

Find  u*  .  h\  .) 


mininuzino 


J.  lu) 


u  ) dxdy 


(5.1) 


among  all  functions  satisfying 


u 


1 


s  ,  1*1 , . . . #N 

l 


with  ii'***"v  d  Partltion  of 

For  this  problem  an  appropriate  finite  element  method  uses  "tent"  functions  in  C(c) , 
with  support  of  the  form  (7)  : 


Figure  5.1 


If  we  denote  by  z . z  the  vertices  of  the  rectangular  mesh  covering  then  the 

1  n 

element  v^  with  support  center#**!  at  is  determined  by  beinq  piecewise  linear  and  satis¬ 

fying  the  interpolation  conditions: 

v  ( z  )  -  l 

[*>.2)  1 

v.ZiM)  .  v(ziM)  -  V(Z1+KI  -  v(Zi_Kl  -  v(Zi  +  K4l)  =  v  'Zi-K-  j  *  ■  ° 


with  F  the  number  of  points  •  z  n\  each  row  of  the  mesh. 


To  obtain  the  matrix  A  in  (3.5)  we  observe  that  for  the  case  support  (v  t 


(5.3) 


1-) 

I x-3 1  -  i .  !  l-j ;  -  k 

otherwise 


In  case  support  (v  )  f  support  (v  )  n  0  #  ($!  the  entries  of  the  corresponding  :  ■  >s  ; 

calculated  according  to  the  geometry  of  the  domain  . 

The  matrix  E  in  (.1.5)  is  computed  by  considering  the  qeometry  of  the  part  it:.-:. 

In  particular  if  support  (v  )  c.  then  the  elements  of  the  corresi  or..: i r.  • 

IN  1  ) 

E  are  determined  by 


(5.-1) 


If  .-1 
if  ipj 


Since  for  l  such  that  support  (v^)  a  =*  (<•)  there  is  no  corresponding  row  m  A 

the  diagonal  of  A  consists  of  positive  entries  only.  Therefore  the  computation  of  t 
solution  of  the  linear  system  (3.6)  can  be  performed  by  one  of  the  iterative  schemes  a 
in  16).  For  the  computation  of  the  numerical  example  to  be  presented,  we  use  a  genera 
JOR  scheme  obtained  by  s| lifting  the  matrix  of  the  system  (i.6)  in  the  following  wa-  : 


wit):  B  *  -  diag  A,  0  *  H-A,  0  -  .  -  1.  This  splitting  corresponds  to  the  s;  littin-:  , 

u) 

matrix  A  according  to  the  classical  JOR  method  (121.  The  convergence  of  rhi*-  itr: 
scheme  to  the  solution  of  the  system  (!.<■»)  from  any  initial  quess  is  guaranteed  hv  The 
in  (o)  »  since  A  and  F  satisfy  the  following  requirements:  A  is  symmetric  nor, nog  a 
definite  with  positive  entries  along  the  diagonal;  the  only  common  vectoi  to  the  null 
of  A  and  K '  is  x  =  d. 


The  computation  of  each  iterant  is  accomplished  in  three  steps: 


T 


( k) 

(1)  Given  b  compute  x  from 

-1  .  (k) 
x  =  B  Cb 

(2)  Solve  for  the  linear  system  of  order  N(<<n) 


(ElB_1E)Y(k+1) 


(3)  Compute  b*k+*'  from 


E'x  -  s 


(k+1) 


-1  (k+1) 

x  -  B  Ey 


Hence-  the  computation  of  b*k+^'  involves  two  multiplications  of  a  vector  by  the  natr.x 
B  *.  This  computation  requires  2n  multiplications  since  B  is  diagonal.  Moreover,  the 
matrix  E'B  ^E  of  the  order  N  <<n,  can  be  initially  factorized  into  the  product  LI ‘ , 
where  L  is  lower  triangular. 

Then  each  step  (2)  involves  only  one  forward  and  one  backward  substitution,  each  of  order 
N.  For  the  actual  computation  it  is  not  necessary  to  store  the  matrices  B,  C,  E ,  but  t. 
store  sufficient  information  for  the  performance  of  multiplication  of  each  of  these  natric.  • 
by  a  vector. 

Under  the  restrictive  assumption,  that  all  boundaries  of  lie  alone  i-.cri- 

1  N 

zontal  and  vertical  mesh  lines,  the  determination  of  the  elements  of  B,  C  and  E  djriru:  t:>. 
computation  is  considerably  simplified.  It  is  sufficient  to  store  an  array  of  dimension  n  ■ 
with  each  row  indicating  the  subregions  which  contain  the  four  parts  of  the  sui  port  of  t h- 

corresponding  finite  element  (these  piarts  are  denoted  by  1,  2,  3,  4  in  Fuiure  5.1).  Ti.,s 

information  contains  also  the  geometry  of  the  domain  .  ,  if  an  additional  romon  ,  . 

assumed  to  contain  all  those  parts  of  the  n  supports  which  are  not  contained  in  ^,...  ... 

The  convergence  of  this  finite  element  scheme  as  h  »  0  to  the  solution  of  ('-.11  •'  i  .  . 

from  Theorem  4.1,  in  view  of  the  convergence  properties  of  this  scheme  when  a;  i  lie.)  t  ■  : 

value  problems  of  the  form  |7|  : 

Aj(u,v)  =  J  fv  dxdy  for  any  v  !!*(  ) 

(5.5) 

/  v  =  0 
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*  u  .» a 

0  to  the  solution  of  (5.51  is  0(h),  if  — ,  7 —  are  hipschitz 

lx  'V 


witn  f  L~  i  <  i  -  0.  It  is  snowr.  in  |7)  that  the  L*'!.)  convergence  of  this  finite 
element  seneme  as 
cunt 1  mtous . 

We  com;  lete  this  section  by  [resent l ng  a  numerical  result  whicn  demonstrates  the  conver¬ 
gence  of  tins  finite  element  scheme,  as  h  -  0,  to  the  solution  of  a  problem  of  the  form 
(5.1).  Trie  [  roblem  we  will  consider  has 


* 

[0,1]  »  [0,11 

1 

E j. 1 1  *  E-,11 

2 

10, h  »  [0,11 

3 

i4.n  «  io, h 

t.  *- 

S3  *  2  • 

In  Table  5.1  we  present  the  values  of  vh  and  an  estimate  of 


h 

v  -  u* 


at  the  points 


lx,y)  »  (4,0)  and  (0,4).  Table  5.1  suagests  that  for  this  problem  the  convergence  of  the 

4.  4. 


finite  element  scheme  is  at  least  0(h) ,  and  possibly  0(h  ) ,  as  h 


0. 


h'1 

T - 

!  ( x  ,  y )  - 

]  h 

..  _V 

■vh-u*'  i 

( x  ,  y )  * 

h 

V 

'k”.. 

v  -  u  *  ! 

8 

I  9.540  | 

.232  | 

17.039 

.040 

16 

!».m  | 

.  06 1  1 

l-’.oio 

.011 

2J 

!  9.747  ; 

.025  j 

17.004 

.005 

32 

!  9.761 

.011  | 

17.002 

.003 

Tahir  ri .  1 
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